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DISJOINT SPARSITY FOR SIGNAL SEPARATION AND 
APPLICATIONS TO HYBRID INVERSE PROBLEMS IN 
MEDICAL IMAGING 


GIOVANNI S. ALBERTI AND HABIB AMMARI 

Abstract. The main focus of this work is the reconstruction of the signals / 
and gi, i = 1,..., N, from the knowledge of their sums hi = f + gi, under the 
assumption that / and the gi s can be sparsely represented with respect to two 
different dictionaries Af and A g . This generalizes the well-known “morpholog¬ 
ical component analysis” to a multi-measurement setting. The main result of 
the paper states that / and the gi s can be uniquely and stably reconstructed 
by finding sparse representations of hi for every i with respect to the concate¬ 
nated dictionary [Af,A g \, provided that enough incoherent measurements gi 
are available. The incoherence is measured in terms of their mutual disjoint 
sparsity. 

This method finds applications in the reconstruction procedures of several 
hybrid imaging inverse problems, where internal data are measured. These 
measurements usually consist of the main unknown multiplied by other un¬ 
known quantities, and so the disjoint sparsity approach can be directly ap¬ 
plied. As an example, we show how to apply the method to the reconstruction 
in quantitative photoacoustic tomography, also in the case when the Griineisen 
parameter, the optical absorption and the diffusion coefficient are all unknown. 


1. Introduction 

Hybrid, or coupled-physics, inverse problems have been extensively studied over 
the last years, both from the mathematical and the experimental points of view. A 
hybrid imaging modality consists in the combination of two types of techniques, one 
exhibiting the high contrast of tissues and a second one providing high resolution. 
Thus, the main drawbacks of the standard imaging modalities can be overcome, at 
least theoretically. Many combinations have been considered, such as optical and 
acoustic waves (photoacoustic tomography [54 ), electric currents and ultrasounds 
(ultrasound modulated EIT 0) or microwaves and ultrasounds (thermoacoustic 
tomography m)- The reader is referred to d isa nm is im E3 for a review on the 
mathematical aspects related to hybrid imaging problems. 

In general, the inversion for such problems involves two steps. In the first step, 
an inverse problem related to the high resolution wave provides certain internal 
measurements. Such internal data are usually functionals of the unknown param¬ 
eters and of certain solutions of partial differential equations (the unknowns are 
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normally the coefficients of the PDE). In the second step, the unknown parameter 
has to be reconstructed from the knowledge of the internal measurements. This 
is sometimes referred to as the quantitative step, since the information on the tis¬ 
sue properties contained in the internal data is only qualitative. In this paper we 
suppose that the first inversion has been performed, and focus only on the second 
step. 

The quantitative step is normally solved with PDE-based methods, by combining 
the internal data with the PDE modeling the problem. Such approach is sometimes 
very powerful in the reconstruction ED [13 ED3 El 0. However, there may be diffi¬ 
culties in using these methods. First, the PDE model may be accurate only in some 
circumstances but not in others m- Second, even if the PDE model is accurate, 
there may be too many unknowns to have unique reconstruction m • Third, even 
in cases when the reconstruction is unique, this may require the differentiation of 
the data P3HHEE], which is known to be an unstable process, or may require 
additional assumptions to be satisfied mmmm®- 

The main focus of this paper is an alternative approach to such problem based 
on the use of sparse representations, as it was first done by Rosenthal et al. in 
quantitative photoacoustic tomography (QPAT) |64]. The internal data in a domain 
D can be often expressed as the product of the unknown(s) and an expression 
involving the solutions of the PDE. (For example, in QPAT the internal data have 
the form H = T/iiq where T is the Griineisen parameter, is the optical absorption 
and u is the light intensity.) Taking the logarithm, the inversion corresponds to 
recovering two functions / and g from the knowledge of their sum 

h(x) = f(x ) + g(x), xgO. 

This problem is, in general, clearly unsolvable. However, it is possible to exploit 
the different levels of smoothness of / and g. Indeed, since / represents a property 
of the medium, such as the log conductivity, it is typically highly discontinuous. 
On the other hand, g is an expression involving the solutions of a PDE, and as 
such enjoys higher regularity properties. As a consequence, / and g have different 
features, and this can be used to separate them by using a sparse representation 
approach. 

Two signals /, g G M n can be reconstructed from the knowledge of their sum 
h = / + g provided that they have different characteristics. More precisely, they 
need to be sparsely represented, i.e., with few atoms, with respect to two incoherent 
dictionaries Af and A g . This method is usually called “morphological component 
analysis” (MCA), and was introduced by Starck et al. in [HH] (see [25l l45l l44l l50l 
for related works and ESJ for a survey on the topic. In 
particular, Donoho and Kutyniok HU 123 first provided a theoretical foundation of 
geometric image separation into point and curve singularities by using tools from 
sparsity methodologies. If compared to these works, the novelty of this paper lies 
in the particular structure of the measurements, as we now describe. 

In this work, motivated by hybrid imaging techniques, where multiple measure¬ 
ments with the parameters fixed can be taken, we extend this method to a multi¬ 
measurement setting. In general terms, this corresponds to the reconstruction of / 
(and gi) from the knowledge of their sums 


hi = f + gi 


i = 1,... ,N. 
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We prove that the MCA approach gives unique and stable reconstruction, provided 
that enough incoherent measurements gi s are available. The incoherence is mea¬ 
sured in terms of their mutual disjoint sparsity. In vague terms, the atoms from A g 
used to represent gi should change for different measurements (see Definition [l] for 
the precise conditions). Numerical simulations show that taking several solutions 
to the relevant PDE yields the necessary incoherence. 

As an example, we discuss the inversion for QPAT, both in the simpler case when 
T = 1 and in the case with non-constant T in the diffusive regime for light propa¬ 
gation. For the T = 1 case, this method has the advantages of being very robust to 
noise and of not requiring a particular model for light propagation, if compared to 
the PDE-based approaches EUHE]. In the case when T ^ 1 there is no uniqueness 
in the reconstruction, even with multiple measurements mb if the parameters are 
piecewise constant, uniqueness can be guaranteed, but the inversion may be very 
sensible to noise m- We propose a combination of the disjoint sparsity signal 
separation method and of the PDE method, which provides a satisfactory recon¬ 
struction, without requiring piecewise constant parameters. Numerical simulations 
are provided. 

This work is structured as follows. In Section |2] we recall the basic notions re¬ 
lated to sparse representations and present the method of morphological component 
analysis. In Section [3J the signal separation method based on multiple measure¬ 
ments and disjoint sparsity is described in detail and the main reconstruction result 
is proved. The numerical implementation of the method, together with the possible 
choices for the dictionaries Af and A g in hybrid imaging, are discussed in Section [4] 
Next, this method is applied to hybrid imaging in Section [5j and several numerical 
simulations are provided. Some concluding remarks are contained in Section [6j The 
proofs of the uncertainty principles stated in Section [3] are given in Appendix [A| 

2. Sparse representations and morphological component analysis 

In this section, we introduce the basic notions related to sparse representations 
and morphological components analysis. 

2.1. Introduction to sparse representations. This presentation follows m • 
Let n be the dimension of the signal space, namely we shall consider signals / £ M n . 
Since in this paper we shall consider only images, we should think of n as being the 
resolution of the image, namely n = d x d, where d is the number of pixels in each 
row and column. However, in this section we shall use the more general notation 
/ £ M n , and think of a signal as a column vector of length n. We now discuss how 
a signal can be represented as a superposition of given atoms in a fixed dictionary. 
More precisely, let A £ M nxm be a dictionary, namely a collection of m atoms, that 
are the column vectors of A. We assume m>n and that A has full rank. Thus, it 
is always possible to express / as a linear combination of these atoms, i.e. to write 

(1) f = Ay 

for some vector of coefficients y £ M m . The most common situation is when m = n 
and A is an orthonormal basis: in this case the coefficient y is uniquely determined. 
However, the situation we are interested in is when m > n. In this case the 
dictionary A is redundant, since / can be represented in many different ways as a 
combination of the atoms in A. In other words, the system 0 is underdetermined 
and has many solutions y £ M m . 
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The key observation is that this non-uniqueness can be exploited by selecting 
the best representation y. One way to measure the quality of a representation y is 
given by its sparsity, which can be quantified by the number of non-zero coefficients 
of y 

IMIo := #( a € {1,..., m) : y(a) ^ 0), 

where the symbol # denotes the cardinality of a set. The representation 0 is called 
sparse if ||y|| 0 <C m. From the theoretical point of view, the sparsest representation 
can be found by minimizing the following problem 

(2) min ||y|| 0 subject to Ay = f. 

yeM m 

The practical search for the minimum poses highly non trivial difficulties, and the 
description of the main issues goes beyond the scope of this work. Algorithms such 
as Matching Pursuit m or Basis Pursuit [36] can often be successfully used to find 
the sparsest solution. More details will be given in Section [4] 

The choice of the dictionary A clearly plays a fundamental role in this context. 
Indeed, a signal / admits a sparse representation with respect to a dictionary A if / 
can be written as combination of few atoms in A. Therefore, the dictionary A has 
to be chosen to capture the main features of the signals we consider. Many choices 
of dictionaries for images are available, and a detailed discussion is presented in 
Section [4j 

In the presence of noise, it is helpful to consider a relaxation of § and to allow 
a small error between the signal / and its representation in terms of the atoms of 
A. Thus, the minimization problem becomes 

min ||y|| 0 subject to ||Ay - f || 2 < e, 

for some 5 > 0, or equivalently 

min ||y|| 0 + A \\Ay — f\\ 2 , 
yeR™ 

for a suitable Lagrange multiplier A > 0. 

2.2. Introduction to morphological component analysis (MCA). One of 

the relevant applications of sparse representations is related to the separation of a 
signal into its constitutive components, provided they have different features. We 
shall describe the method discussed in [69]. Suppose that a signal h E M n can be 
written as a sum 

h = f + g, 

with f,g E W 1 . The problem under consideration is the reconstruction of / and g 
from the knowledge of h, under the assumption that / and g have distinctive char¬ 
acteristics. This assumption can be expressed in terms of sparse representations. 
Namely, suppose that there exist two dictionaries Af E M nxm/ and A g E W ixrn 9 
such that: 

(1) / can be sparsely represented with respect to Af but not with respect to A g \ 

(2) and g can be sparsely represented with respect to A g but not with respect to 

A f- 

Under these heuristic conditions (which will be made precise below), a strategy to 
find / and g may be to find a sparse representation y — \y f g ] of h with respect 
to the concatenated dictionary A = [Af,A g \ and then to write / = Afyf and 
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g = A g y g . As we have seen before, the sparse representation y is the minimum of 
the minimization problem 

min ||y|| 0 subject to [A f A g ] [^] = h, 

yeR m f +m g 

or, in the presence of noise, of 

min ||y|| 0 subject to ||[^/ A g ] [ v v f g ] - h\\ < e. 

Even though this procedure is successful in many practical cases [69l l48l 164 , a 
proof of the correct reconstruction, i.e. / = / and g = g, is only valid in an ideal 
situation, which we now describe. In the case when Af and A g are both orthonormal 
sets, the proof is based on the following uncertainty principle ga Theorem 1]. 

Proposition 1. Let {ai,..., a mA } and {6i,..., b mB } be two orthonormal sets of 
vectors in W 1 . Take a vector h G M n \ {0} and suppose it has the following repre¬ 
sentations 

h = Ay A = By B 

with respect to A = [ai,..., a mA ] and B = [b\, . . . ,b mB ]. Then 

WvaWo + Wvb 11o > 2 /-V 

where 

(3) M = max |(a*, bj) 2 \ 

hi 

is the mutual coherence. 


We provide a proof for completeness. 

Proof Let Sa and Sb denote the supports of yA and yB , respectively. By assump¬ 
tion we have ys(j) = (h, bj) 2 = ^Z ie s A VA{i){ai,bj) 2l whence by Cauchy-Schwartz 
inequality we obtain |y B (j)| < \\va \\ 2 VWSaM. Thus, since ||?m|Io = #S A and 
\\Vb\\o = # 5 'b, this implies 

l|ys|| 2 < \\va \\ 2 V(#Sa)(#S b )M < H2/AH2 (Iloilo + Ill/sUy- 

This concludes the proof, since the orthonormality of A and B implies ||^a|| 2 = 
||ft|| 2 = \\Vb\\ 2 (Lemma |T] part 1). □ 

In general, it is easy to see that if A and B are orthonormal bases then 1 />Jn < 
M < 1 [45]. As a simple consequence of this result ga Theorem 2], we have 
that if y 1 G M 2n and y 2 G M 2n are two representations of h with respect to the 
concatenated dictionary A = [Af,A g \, then 

||y 1 || 0 + ||y 2 || 0 >2/M. 

Therefore, if / and g have representations yf and y g satisfying ||^/|| 0 + ll%llo ^ , 

then the above method provides the correct reconstruction. 

In practice, the assumption \\yf\\ 0 + ||%|| 0 < 1/Af is almost never satisfied, and 
so the above argument remains only a theoretical speculation. However, when 
the multi-measurement case is considered, correct and stable reconstruction can 
be rigorously proved. This theoretical result is also validated by several numerical 
simulations. These aspects are discussed in the following sections. 
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3. Disjoint sparsity for morphological component analysis 

3.1. Introduction and main assumptions. Motivated by several hybrid imag¬ 
ing modalities (see Section^, we generalize the MCA problem to a multi-measurement 
setting. The reader is referred to da [sni mi for other similar variations. 

Let /o,..., h]\f G M n be TV signals that can be decomposed as 

hi = f + <ji, i = 1,... ,7V, 

with f,cji G M n . We want to study the problem of finding / and gi from the 
knowledge of hi for i = 1,..., TV. The case TV = 1 was discussed in the previous 
section. We shall show that as TV becomes bigger, the above problem becomes much 
more treatable, and that the sparsity approach introduced before always provides 
the correct reconstruction, also in the presence of noise. As before, let Af and 
A g be the dictionaries with respect to which / and gi have sparse representations, 
respectively. Note that all the gi s can be sparsely represented with the same 
dictionary A g : this is a crucial assumption of this approach. Assume that the 
atoms of Af are normalized to 1 and that 

(4) the atoms of A g constitute an orthonormal set of M n . 

Thus, A g can be completed to an orthonormal basis [A g ,A^] for some Aj G 

The reconstruction method applied to this case consists in the minimization of 

(5) min ||j/|| 0 subject to \\A f y f + A g y l - h t II < e, i ** 1 N, 

ye m.™f+ N ™9 y 

where we have used the notation y = 'y^]- Here, the superscript t 

denotes the transpose. To model the case with added noise, we write 

(6) hi = f + gi+ni, 

where / and the gi s represent the true signals, the hi s are the measured signals 
and rii is such that 

(7) ||rai|| 2 <?7, i = 1 , ... ,N 
for some small 77 > 0. 

In the applications we have in mind (Section [ 5 ]) , the signal / represents (the 
logarithm of) a physical constitutive parameter, while the g^s usually quantify the 
injected fields, e.g., the electric field or the light intensity. As such, / is given and 
fixed, and we have no control on it. On the other hand, the g^s come from the 
measurements, and can be indirectly controlled. More precisely, the g^s depend on 
the solutions of a certain PDE, whose coefficients are unknown, but whose boundary 
values can be chosen: in this sense the gi s can be controlled. It is therefore natural 
to give some assumptions on the gi s. 

The main requirements are that the gi s should be sufficiently many and incoher¬ 
ent This will be mainly expressed by means of their disjoint sparsity with respect 
to A g . (Disjoint sparsity was used in [24], while joint sparsity has been extensively 
used in compressive sensing US H>].) We shall therefore write 

(8) —/|| 2 < pj, || A g y g — gi || 2 < Pgi i = l,...,N 


1 Similar assumptions of enough indepen dent measurements are required also when using PDE 
methods for hybrid inverse problems (see § 


5.1). 
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for some pf,p g > 0. The approximation allows a small error between the true signals 
and their sparse representations. We shall prove that under suitable assumptions, 
a minimizer of © provides the correct reconstruction, up to a factor that is small 
in 5 := pf + p g + y. In terms of the coefficient vectors {; y g ,..., y^}, the required 
assumptions can be written as follows. 

Definition 1. Take /3,D > 0, N G N*, yf G M m/ and G M m ^. We 

say that {y g ,..., y ^} is a (/3, D)-complete set of measurements if the following two 
conditions hold true: 

CS1: if | yg(a) — y 3 g (a)\ < /3 and y g (a)y J g (a) ^ 0 for some a G {1,... ,m g } then 

i = j m , 

CS2: for every q G M m/ such that ||^4/^|| 2 > D and ||M.j-A/g|| 2 < 2/3 there holds 

N N 

(9) (supp Q \ SU PP Vf) + 'y ] 7 \ supp y g ) > # U suppi/* + H^llo , 

i=1 i= 1 

where S q = {a : \( t A g Afq)(a)\ > 1}. 

Let us comment on these two requirements. The first condition, CS1, is a con¬ 
straint on the incoherence of the <^’s in terms of the values of their coefficient 
vectors. More precisely, the coefficients relative to the same atom for two different 
Pi s cannot be too close. This assumption could be relaxed by allowing a fixed 
number of y g to have the same coefficients, but for simplicity of exposition we have 
decided to omit this generalization. 

While CS1 is mainly a technical hypothesis, condition CS2 is the main assump¬ 
tion related to the disjoint sparsity of the cji s. Indeed, when the sets supp y g are 
small (the representation is sparse) and substantially change when i varies (dis¬ 
joint), then it becomes possible to satisfy CS2 by taking N large enough, since the 
right-hand side is bounded by \\yf\\ 0 + m g . Note that CS2 can always be satisfied 
by choosing the y g so that #{i : y l g (ot) = 0} > #(Ji su PP^ + \\Vf\\o f° r every cq 
but this represents only a worst case scenario. In general, the smaller supp yf and 
suppi/* are, the easier it becomes to satisfy CS2. 

Remark 1. It is expected that the number of measurements N should increase as 
the noise level y becomes bigger. Indeed, to have a good reconstruction, we need 
to keep 5 = pf + p g + y small. Thus, if the noise level y increases, the quantities 
pf and p g have to become smaller. In other words, the sparse approximations in 
© have to be more precise, which in turn requires the support of yf and y g to be 
bigger, and so a higher number of measurements is needed to satisfy CS2, as we 
heuristically observed above (see Example [l] below for a concrete example). 

3.2. Uncertainty principles. As we have already pointed out, the smaller supp yf 
and supp y g are, the easier it becomes to satisfy CS2. Moreover, when Af and A g 
are two orthonormal bases, by Proposition [l] we have 

II 4o + II^A/^Ho ^ 2 /^ 

where M is the mutual coherence of Af and A g defined in <©• However, this 
inequality is not directly applicable to our context, since 11^4/11 0 does not appear 
explicitly in (|9| and A g will not be a basis in the applications. The following 
generalization of the uncertainty principle guarantees that the same bound holds, 
provided that D is big enough. The proof is given in Appendix | A. 1| 
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Proposition 2. Assume that Af is an orthonormal basis of W 1 and that A g is 
an orthonormal set ofW 1 and let M denote their mutual coherence. There exists 
D > 0 such that 

H<4 + #{<* : | (*A g A f q)(a)\ > 1} > 2/M, 
for every q E M n with ||Ayg|| 2 > D and ||M.j-A/<z|| 2 < 2/3. 

Remark 2. In the single measurement case (if Af is an orthonormal basis), CS1 
and CS2 correspond to the condition \\ijf\\ 0 + ||%|| 0 < 1/M, thereby reobtaining 
the hypothesis discussed in the previous section. Indeed, if N = 1 then CS1 is 
immediately satisfied. Moreover, CS2 becomes 

#supp<?\suppy / + #{a: \{ t A g A f q){a)\ > l}\suppy s > \\y g \\ 0 + ||j//|| 0 . 

By Proposition |2j the left hand side is bounded by below by 2/M — (||j//|| 0 + ||y 9 || 0 ), 
and so the above inequality is satisfied provided that ||^/|| 0 + ll%llo < 

The more incoherent the two bases are, the bigger 2/M is, and therefore the 
bigger the left-hand side in (|9| is. It is not a surprise that the incoherence of the 
bases makes these assumptions easier to be satisfied, since this was the starting 
point of the MCA discussed in the previous section. 

Let us discuss a simple example to show the relation between the sparsity of the 
coefficients and the number of measurements N. 


Example 1. Let us consider the simplest problem of the separation of spikes and 
sinusoids in ID. Let Af = I n E M nxn be the identity matrix and A g be the Fourier 
basis. Their mutual coherence is M = 1 / y/n. Let yf E M n be such that \\yf \\ 0 = k 
and y g E M n be such that ||^|| 0 = l and with disjoint support, for i = 1,... ,7V. 
We would like to investigate the link between the number of measurements N and 
the sparsity of yf and y g , i.e. with the quantities k and /. 

Thanks to the assumption on the supports of the y g ' s, condition CS1 is automat¬ 
ically satisfied for any /3 > 0. As far as CS2 is concerned, unfortunately we cannot 
check its validity for all possible choices of q. Thus, we make an heuristic choice 
with one of the worst possibilities, the Dirac comb, for which the inequality of the 
uncertainty principle becomes an equality. Set qs(od) = 1 when a is a multiple of y/n 
and qs{a) = 0 otherwise. It turns out that A g qs = and so ||^|| 0 = #S qs = ^Jn. 
Provided that l < y/n, it turns out that condition CS2 for y$ is satisfied if 


2 k 

tvTT 


T / < \fn- 


Thus, the bigger the number of measurements is, the bigger k and l can be. 

Assume now that A g contains only a subset of the Fourier basis of cardinality 
m g (as we shall do in the numerical simulations), and select a vector q$ such that 
Mo = = yjn as before. Thus, CS2 is satisfied if 


2k < N(y/n — l) + y/n — m g , 


independently of the supports of the y l g s. For a fixed l < y/n, the higher N is, the 
bigger k can be. 


When considering wavelets and the Fourier basis, as in the numerical simulations 
of this paper, there is the added difficulty that the mutual coherence is not small, 
but in fact of order 1, which makes the uncertainty principle discussed above of no 



DISJOINT SPARSITY FOR HYBRID IMAGING 


9 


practical use. However, if we consider Haar wavelets and low frequency trigono¬ 
metric polynomials, as in Section [5j it is possible to improve this bound (a similar 
result for spikes and sinusoids is given in m Lemma 1]). We now discuss this 
result. 

We consider two dimensional signals in M n , where n = dxd and d = 2 J for some 
J G N*. Let Af be the associated orthonormal basis consisting of 2D periodized 
Haar wavelets (see § A.2). Let A g consist of low frequency non-constant trigono¬ 
metric polynomials. More precisely, for L = 1,..., d/2 — 1 consider the following 
families of real sinusoids 

Xi( a ) = c i sin(27r/iaq/d) sin(27r/2<T2/d), h,h = 1,..., L, 

xf( a ) = c 2 sin(27rZiQ'i/d) cos(27r/2<T2/d), l\ = 1,..., L, I 2 = 0,..., L, 

xf( a ) = c f cos(27r/iGi/d) sin(27r/2<T2/d), l\ = 0,..., L, I 2 = 1,..., L, 

Xi (a) = cf cos(27rZiGi/d) cos(27r/ 2 G2/d), l G {0,..., L} 2 \ {(0,0)}, 


( 10 ) 


where cj > 0 are suitable normalization factors chosen so that ||xj|| 2 = 1. Let A g be 
the collection of all these real sinusoids. Assumption @ is satisfied. The number 
of atoms is given by m g = 4L 2 + 4 L. The proof of the following result is given in 
Appendix |A.2| it is heavily based on the structure of Haar wavelets, and so it is 
expected that this result would fail with smoother wavelets. 


Proposition 3. Let Af and A g be the dictionaries of 2D Haar wavelets and low 
frequency non-constant real sinusoids discussed above, respectively . Assume that 
L < 2 b for some B < J — 2. There exists D > 0 such that 

Ikllo > E (2 J - j -2 L)\ 

3 =1 

for every q G M n with ||Ajg|| 2 > D and || t A^Afq\< 1 < 2/3. 

Example 2. In Section [5] we shall set J = 7 and L = 15. In this case, the 
above inequality reads ||g|| 0 > 1160, which is sensibly better than the uncertainty 
principle based on the mutual coherence. Thus, arguing as in Remark [2j in the 
single-measurement case, condition CS2 is satisfied provided that 

(11) 2||^|| o + ||%|| o <1160. 

It should be observed that in the numerical simulations we add also the constant 
matrix to the dictionary A g . Even if this is not allowed by the above result (only 
4 wavelets are needed to represent it), it seems not to create any issues for the 
reconstruction. This might be due to the following remark: most images are not 
constant. 


Remark 3. We conclude this part by observing that these uncertainty principles 
always take into account the worst case scenarios; namely, for most vectors the 
minimum is not attained. Improved bounds (called robust uncertainty principles) 
that hold for the overlwhelming mojority of vectors were proved in |32| for spikes and 
sinusoids. Moreover, it was proven that the separation problem can be succesfully 
solved in most cases by using Z 0 minimization, provided that corresponding sparsity 
conditions are satisfied. These conditions are much weaker than the ones based on 
exact uncertainty principles. It would be interesting to investigate whether such 
extensions hold in our context too. 
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In view of this, while the uncertainty principle in Proposition [3] does not contain 
the term corresponding to S q , it seems reasonable to assume that for most vectors 
the quantity #S q is not small when ||g|| 0 is close to the minimum. This explains 
why CS2 is easily satisfied when N is bigger, at least for most q. 

3.3. Main result. The main result of this section states that if the y l g s constitute 
a complete set, then the signals / and cji can be stably recovered from the knowledge 
of hi m f + rii by minimizing 0- 

Theorem 1. Assume that 0 holds true. Let /?, D, pf, p g ,g > 0 and N G N* be 
such that e := pf + p g + rj < /3/3. Take yf G M m/ and let {y g , ..., y^} be (/?, D)- 
complete. Assume that f\gi,ni,hi G M n satisfy 0, 0 and (|8| and let yf G M m/ 
and y g G M 7 ™ 9 realize the minimum of 


(12) min \\y\\ 0 subject to || A f y f + A g y % - hi\\ < e, i = 1. 

yeR m f +9 

Then for every i = 1,..., TV 

\\AfVf ~ f || 2 < (3^ + l)e> ||Agyg - gi\\ 2 < (3 D + 2)e. 




Thanks to this result, the reconstruction can be performed by minimizing (12) 
and then taking f & Afyf and gi « A g y g . The practical details of such minimiza¬ 
tion will be discussed in Section 01 

In the proof of this theorem we shall make use of the following properties, whose 
proofs are immediate. 


Lemma 1. Let A G M nxm constitute an orthonormal set of M n and let A 1 - G 
^nx(n-m) com pi e f e a to an orthobasis o/M n . The following properties hold true. 

(1) For every v G M m we have ||Ai;|| 2 = |H| 2 - 

(2) For every v G M n we have ||Mp|| 2 < ||p|| 2 . 

(3) We have l AA — I and l A^A = 0. 

(4) For every a,b G M m we have 

||a + 6|| 0 = IMI 0 + #(supp b \ supp a) - #{a : a(a) = -b(a) / 0}. 


We are now ready to prove Theorem [1] 

Proof of Theorem [7} Write / := A f y f , g { := A g y l g , 

:= — l A g Afef and r 2 := e g — e g . Since and y g satisfy the constraint in (12) 
we have 


e f : = Vf ~ vp e g : = y g - y g 


II A fVf - A fVf + A a y l g - Agtgh < II ( A fVf + A gV l g - h i) + {K - A g y l g - A f y f )|| 2 

< 6 + 11 / - A fVf II2 + II 9 i - A g y g II2 + 11^2 11 2 

< 2e 


where the last inequality follows from 0 r and |8|. As a consequence, since r l = 
t A g (Afyf — Afyf + A g y g — A g y g ) (Lemma [l] part 3), by Lemma |T] part 2 we obtain 
||r*|| 2 < 2e, whence 

(13) IkiL < 2e. 

Moreover, by Lemma [l] parts 2 and 3 we obtain 

(14) llTj- A e /ll2 = - A fVf + - A gV l g) II 2 < 2e - 
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Since yf and y l realize the minimum of ( 12 ) we have ||?/|| 0 < ||f /|| 0 or, equivalently, 


N 


(Wvfh - WiffWo) + Hdbsllo - Iloilo) ^ °- 


i=1 


where we have set y = t [ t yf , t yg, ..., t yg]- Thus, since yf=yf + ef and y 1 = 
y* + e*, Lemma [l] part 4 yields 

(15) #(suppe/ \supp yf) - #{a : y f (a) = -e f (a) 7 ^ 0} 

N N 

+ 5Z#( SUppe s \ SU PP^) ~^Z#{ a : y l g( a ) = ~ el g ( a ) 7 ^ 0} < °- 


(16) 


i=1 i=1 

Observe now that by ( |13| ) we have 

#{a : y f (a ) = -e f (a) ± 0 } < ||y /|| 0 , 
supp e g = supp (e g + r l ) D {a : |e fl (a)| > 3e}. 


Since 3e < /? and { fig ,..., y g } is (( 3 . D)-complete, by <E§ and Definition |lj (condi¬ 
tion CS 1 ), we have 
(17) 

N N N 

: yl( a ) = - e l( a ) Y 0} = # (J{a : y l g (a) = -e*(a) ^ 0} < # |J suppy*. 
2=1 2=1 2=1 

Inserting (16) and into p5) gives 

AT AT 

#(suppe / \suppy / )+y]#({a: |(*A g A / e / )(a)| >3e}\suppy*) < ||y/|| 0 +# (Jsuppy*. 




2=1 


Set q = eg/(3s). Since suppej = suppg we have 


N 


N 


#(supp q \ supp yf) + #(S q \ supp y g ) < ||y /|| 0 + # |J suppy*. 


2=1 


2=1 
it j 


where S q = {a : \( t A g Afq)(a)\ > 1}. Moreover, by (14) we have ||Mj-A/-g || 2 < 2/3. 
Thus, since is (/?, D)-complete, by Definition [l] (condition CS2) we 

obtain that /1 | 2 < 3 De. As a result, by construction of ej, (| 8 | and the triangle 
inequality we have 


(18) 


| A f y f - f || < 3 De + Pf < (3 D + l)s. 


This proves the desired bound for /. It remains to prove the corresponding estimate 
for gi = A g y l g . In order to do this, observe that by ( p~ 8 | ) and the triangle 
inequality we obtain 


9i ~ 9i\\ 2 < ||(/ + 9i ~ hi) + (/ — /) + ^i || 2 ^ 6 + + P/ + 9 < (3D + 2)s. 


This concludes the proof. 


□ 
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4. Numerical implementation 


4.1. Orthogonal Matching Pursuit. The simplest available algorithm for the 
minimization of problems of the type 


min llvL 


subject to \\Ay - f || 2 < 5, 


for / G M n and A G M nxm , is the Orthogonal Matching Pursuit (OMP) [59, 27]. 
This algorithm belongs to a wider class of methods called greedy algorithms, in 
which the set of the used atoms of A is increased step by step. In OMP, the best 
coefficients for the atoms are recomputed at each iteration, which makes it more 
efficient compared to the standard matching pursuit. Greedy algorithms, and hence 
OMP, are not a priori equivalent to the minimization of the above problem, and the 
convergence to a minimizer is not always guaranteed, but they have been proved 
to perform well in most cases m- The reader is referred to m for more details 
on this topic. 

The adaptation of OMP to the problem of our interest is quite straightforward. 
By Theorem [l] we need to minimize i.e. 

min llvllo subject to \\A f y f + A g y l g — h t || < e, i=l,...,N, 

y eR rn f + 9 

given N signals hi G M n and two dictionaries Af G W ixrn f and A g G M nxm ^. 
Setting 



A f 

A 9 

0 

... o " 


~Vf~ 

1 


"hi" 

A = 

A f 

0 

A g 


, y = 

yg 

y 9 

and h = 

h 2 





0 






Af 

0 


0 Ag_ 


yf. 


h N _ 


the above problem is equivalent to 

min ||i/|| 0 subject to \\Ay — h \\< VNs, 

y£R rn f + Nrn 9 


for which the OMP can be applied directly, as discussed above. As before, OMP is 
not guaranteed to converge to a minimizer of this functional. In other words, even 
though Theorem [l] gives (almost) exact recovery via the minimization of (12), OMP 
may not provide a faithful reconstruction. However, the numerical simulations 
(Section [4| suggest that in practice a correct reconstruction is always achieved via 
OMP. It would be interesting to consider this issue from the theoretical point of 
view, but this goes beyond the scope of this work. 

Let us also mention other alternatives for the minimization of these problems, 
such as Basis Pursuit [SB], Block-Coordinate-Relaxation [26], Iterative Hard Thresh¬ 
olding [22] and Stagewise OMP [43]. In particular, in Basis Pursuit the ^-penalization 
term is substituted by an i 1 term, which makes the functional convex: the min¬ 
imization can be easily achieved with standard tools. This would complicate the 
corresponding reconstruction result (Theorem [TJ and consequently the assumptions 
on the original signals (Definition [l]) . Achieving this would be a natural follow-up 
of this work. 

The Block-Coordinate-Relaxation method was adapted to the minimization of 
the functional related to the separation problem with N = 1 measurement |69| . 
This method, sometimes referred as MCA, is a combination of Basis Pursuit (the 
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i 1 norm is used instead of the norm) and of a two-step iterative procedure 
where the two components of the signal are minimized separately. Moreover, the 
minimization is expressed directly in terms of the signals Afi/f and A g y g , without 
the need of storing the full dictionary matrices, which may take a lot of memory 
space. It would be very interesting to see whether this method can be generalized 
to the multi-measurement case considered in this paper. 

4.2. Dictionaries for image content. Let us now discuss what choices may be 
done for the dictionaries Af and A g in the context we are interested in, namely 
hybrid imaging inverse problems. As we shall see in Section [5j in such problems 
the signal / will typically represent (the logarithm of) a constitutive parameter 
of the tissue under investigation, like for instance the electric permittivity, electric 
conductivity or the optical absorption. As such, the image / can be supposed to 
be piecewise smooth: the discontinuities are usually the inclusions we would like to 
determine. On the other hand, the signals g^s usually represent the injected fields, 
such as the electric field or the light intensity, and are the solutions to certain 
PDE. As such, they enjoy higher regularity properties, and the images cji s can be 
supposed to be smooth. 

These different features represent the foundation of the signal separation method 
discussed in the previous section. In order to exploit this diversity it is necessary to 
choose suitable dictionaries Af and A g , with respect to which / and the cji s have 
sparse representations, respectively. 

As far as Af is concerned, wavelets have been widely used to sparsely represent 
piecewise smooth images EOJ. In particular, Haar wavelets will be used in this work: 
the choice is motivated by Proposition [3] It is worth noting that in recent years a 
large number of new multi-dimensional transforms have been introduced in order 
to capture the directional features of two-dimensional images [65] . Among the most 
known, there are curvelets m, ridgelets [29] and shearlets [58j|57l[52]. The use of 
these transforms may give better results, but a deep investigation of the best choice 
for the dictionaries goes beyond the scope of this paper. Thus, we have decided to 
make the most classical choice of wavelets, since, even though possibly not optimal, 
it is sufficient to properly illustrate the disjoint sparsity signal separation method. 

The situation is simpler for the choice of A g . Indeed, a good representation of 
smooth functions may be obtained by choosing low frequency trigonometric polyno¬ 
mials. This is a simple consequence of the correspondence between the smoothness 
of a function and the decay of its Fourier transform. This represents our choice for 
A g in this paper. 

These dictionaries purely represent a general guideline for the choices of Af and 
A g . Additional information on the particular physical model may be used to select 
dictionaries more adapted to the features of the images under consideration. 


5. Applications to hybrid inverse problems 

5.1. Introduction. We have seen in the introduction that a hybrid problem usu¬ 
ally involves two steps. First, internal functionals are measured inside the domain 
and, second, from their knowledge the unknown coefficients of the PDE have to be 
reconstructed. These internal data are linear or quadratic functionals of the un¬ 
knowns and of the solutions of the direct problem. Let us mention some examples. 
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• In photoacoustic tomography muni hues] the internal data take the form 

H{pc) = T(x)ii(x)u{x), xGH, 

where T is the Griineisen parameter, g is the optical absorption and u is the 
light intensity. The main unknown of the problem is fi. The photoacoustic 
image H gives only qualitative information on the medium, as fi is multiplied 
by T and u. The problem of quantitative photoacoustics is the reconstruction 
of fi from H. Under the diffusion approximation, the light intensity u solves 

(19) — div(D\7u) + fiu = 0 in U, 

where D is the diffusion parameter. This equation should be augmented with 
suitable boundary conditions, such as of Dirichlet or Robin type. 

• In thermoacoustic tomography m the internal data take the form 

H(x) = a(x)\u(x)\ 2 , xGH, 


where a is the unknown conductivity and u is the electric field. The problem of 
quantitative thermoacoustics is the reconstruction of a from H. In the scalar 
approximation, u is the solution of the Helmholtz equation 

A u + (c d 2 + iuj(j)u = 0 in U, 


where uj is the angular frequency. As before, this equation should be augmented 
with suitable boundary conditions. 

• In microwave imaging by ultrasound deformation Q21E] the internal data take 
the form 

H(x) = £(x)u(x) 2 , xgD, 

where s is the unknown electric permittivity and u is the electric field. The 
problem is now to reconstruct 6 from H. In the scalar approximation, u is the 
solution of the Helmholtz equation 

A u + uj 2 eu = 0 in D, 


augmented with suitable boundary conditions. 

• In ultrasound modulated diffuse optical tomography QH M ES] the internal 
data are di v(r 2 V/j), where u solves (19) and /i is the optical absorption. 


In all the above examples, the measurement H is the product of the desired unknown 
and other quantities. Thus, the problem is extracting the desired unknown from 
H. PDE techniques are usually used to solve the problem, but, as discussed in the 
introduction, have several drawbacks. 

All the above problems consist in the determination of two functions from the 
knowledge of their product. Taking logarithms, in a multi-measurement setting 
this is equivalent to finding f{pc) and gi{x) from the knowledge of their sums 


h(x) = f(x) + gt(x), x e Q. 


The disjoint sparsity signal separation method can be applied since / and the gi s 
have different level of smoothness, and so can be sparsely represented with respect 
to different dictionaries (see § |4.2| ). 

In particular, Theorem |T] guarantees unique and stable reconstruction of the 
the coefficients, provided that we can construct many incoherent rq of the above 
problems (by changing the boundary values), so that the corresponding gi s give a 
complete set of measurements, according to Definition [l] As we shall see below in 
the numerical simulations, the conclusion of Theorem [l] is verified in several cases, 
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choosing different boundary values. Unfortunately, it is not possible to check the 
conditions of Definition [l] numerically, as this would require an infinite number of 
computations. These conditions are used to heuristically inform the choice of in¬ 
coherent illuminations. On the other hand, at the current state, we are unable to 
give a general strategy to build boundary values so that the corresponding solu¬ 
tions will satisfy the hypotheses of Definition []] this represents a very interesting 
open problem regarding the interplay of the PDE theory with the disjoint sparsity 
approach (see Section [6|. 

It is worth mentioning that similar assumptions of incoherence, usually in terms 
of linear independence of gradients of solutions, are often necessary for the PDE- 
based reconstruction methods (see, e.g., [33 mi ini m HI H3 m El ED and references 
therein). 

As an example, in the rest of this section we shall apply the method to the 
reconstruction in quantitative photoacoustic tomography. All the other modalities 
mentioned above can be treated with minor modifications. 

5.2. Quantitative photoacoustic tomography, the case T — 1. In photoa¬ 
coustic tomography, the object under investigation is illuminated with light radia¬ 
tion, whose absorption causes local heating of the medium. The temperature rise 
creates an expansion of the tissue, thereby producing an acoustic wave, that can be 
measured on the boundary of the domain. The first step of this hybrid modality 
consists in the recovery of the amount of absorbed radiation by inverting the wave 
equation, from the knowledge of the boundary values. This problem has attracted 
much attention in the last years: the reader is referred to m and to the references 
therein for more details. In this paper, we shall suppose that the first step has been 
performed, and that we have access to the amount of absorbed radiation 

H(x) = T(x)n(x)u(x ), xgD, 

where T is the Griineisen parameter, fi is the optical absorption and u is the light 
intensity. The problem of quantitative photoacoustic tomography consists in the 
reconstruction of fi from the knowledge of H. Much research has been done on this 
over the last years, see e.g. [39[ |64| [40j HU CM ESI SSI EM [62] and references therein. 
Sparse representations were first used in m- Recently, one-step methods have 
been introduced to perform both steps at the same time ECO IH|: it would certainly 
be interesting to see whether the signal separation approach can be employed in a 
one-step reconstruction method. 

Light propagation may be modeled by a radiative transport equation or, when 
the scattering coefficient is large and the absorption is small, by its diffusion ap¬ 
proximation ( fl9| ). We consider here the simplest case when T = 1, the general case 
is discussed later in § |5.3| By using multiple measurements, fi can be recovered 
both in the transport regime m and in the diffusive regime EU, under suitable 
assumptions on the solutions. 

We now describe how to apply the disjoint sparsity approach to this problem. If 
compared to the PDE approach, the separation of fi and u does not require the use 
of a PDE, and so no specific model of light propagation has to be assumed for the 
inverse problem. Moreover, no a priori knowledge of relevant boundary conditions 
(which may be unrealistic) is required, in contrast to the PDE approach. These 
aspects should make this approach more feasible. For simplicity, we shall construct 
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( 20 ) 


r - a u, 

\ Ui = 


the solutions iq via ([19]) with D — 1 and Dirichlet boundary values, namely 

—A Ui + /mi = 0 in 
(fi on dQ. 

However, this equation will not be used for the inversion. 

The joint sparsity method will be applied as follows. Let fl denote the true 
absorption. After constructing N solutions fti,... ,rat to the above equation, the 
quantities 

Hi(x ) = fl(x)ui(x ), x £ fi, 

are measured. Taking logarithms and adding white Gaussian noise n$, we obtain 


( 21 ) 


hi = log/2 + log ft* + rii, i = 1,... TV. 


The reconstruction of fl from the knowledge of the ft^’s exactly corresponds to the 
problem discussed in Section [3] The method based on Theorem [l] and whose nu¬ 
merical implementation is described in Section [4] will be used for the reconstruction. 

In all the examples, we shall consider the two-dimensional domain [0,1] 2 with 
128 x 128 pixels. As far as the choice for the dictionaries is concerned, we let Af 
consist of the orthonormal basis of Haar wavelets (as in § A.2) and let A g consist 


of 961 low frequency trigonometric polynomials (as in (10), including the constant 
matrix), periodic over [0, l] 2 . 

The choice for A g forces to choose periodic boundary conditions, and so we set 

Lpx{x) = 1 , 

(f 2 (x) = 1 - sin(27rxi)/4, 

(22) Vz( x ) = 1 - sin(27nr 2 )/4, 

y>±{x) = 1 — cos(27tx 2 )/4, 
ips(x) = 1 — cos(27rxi)/4. 

(For physical constraints, all the quantities have to be positive.) Non-periodic 
choices for the boundary conditions would be allowed with no added difficulties for 
the reconstruction, provided that the dictionary A g is properly chosen. The above 


boundary values are expected to generate incoherent solutions to (20) in such a 
way to satisfy the conditions of complete sets as in Definition [l] In this case, the 
assumptions of Theorem [l] would be verified and the disjoint sparsity separation 
method would be guaranteed to provide the right reconstruction, even in presence of 
noise. However, this is not guaranteed a priori: as mentioned above, the conditions 
of Definition [l] are very hard to check. 

5.2.1. Example 1 - convex inclusions. We consider three convex constant inclusions 
in a homogeneous background, as shown in Figure [la] We choose to stop the 
iteration procedure of OMP after 1500 iterations, which gives satisfactory results. 
More accurate stopping criteria may be considered [64], but this is not among the 
aims of this work. 

In a first experiment we consider one noise-free measurement, namely we set 
N = 1 and n% = 0 in ©. The results are shown in Figure [l] The solution 
to (20) with boundary value pi = 1 is shown in Figure [lb] and the corresponding 


measurement Hi = /Hq in Figure [lc] As it is evident from the images, the inclusions 
are still clearly recognizable in Hi , but the corresponding values of the absorption 
are corrupted by the multiplication by ui. The sparse separation approach is thus 












DISJOINT SPARSITY FOR HYBRID IMAGING 


17 



(a) The true absorption p 



(c) The datum Hi 


Figure 1. Example 1. Noise-: 


0.98 

0.96 

0.94 

0.92 

0.9 

(b) The true intensity u\ 



(d) The reconstructed fi 
case with one measurement. 



applied as discussed above, and the reconstructed p is shown in Figure |ld| The 
relative reconstruction error is 

l|1 ° g ,^~?°f /i|l » * 1.5%. 

l|logWl 2 

This shows that, in absence of noise, the reconstruction is almost exact, even with 
only one measurement. This is in total agreement with the theoretical result dis¬ 
cussed above. In the absence of noise (77 = 0), it is possible to satisfy the conditions 
in Definition [T| in view of Proposition [3] and Example [2] Indeed, the sparse approx¬ 
imations of log p and of log ft 1 with 500 wavelets and 100 sinusoids, respectively, 
give small approximation errors, namely pf « 0.54 and p g ~ 0.05. By jn}, CS2 
is satisfied. Thus, the assumptions of Theorem [l] are verified with e « 0.59, and 
the resulting recontruction is very good (||log/i — log/r|| 2 ~ 0.29). As already men¬ 
tioned in Example [2] this argument is not fully rigorous, since Proposition [3] does 
not allow the constant image to be in the dictionary A g . 

In a second experiment (Figure [2] we add white Gaussian noise n* in ( [TT] ). The 
noise level is so that 

|n ll ?-- 2 -. || « 17.6%. 

11 log WUi )\\ 2 

We tested the reconstruction procedure for TV = 1,..., 5 and the boundary values pi 
as in (22). The data Hi, i = 1, 3, 5, are shown in Figures [2a| [2c] and [2eJ respectively. 
The reconstructed /i’s for N = 1, N =5 3 and N = 5 are shown in Figures [2b] [2d] and 
[2f[ respectively. The reconstruction errors for N = 1,..., 5 are given in Table IT] 
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(a) The datum Hi. 



(c) The datum if 3. 



(e) The datum H$. 



(g) if 2 with TV 



(h) 11 with TV, N = 2. 


Figure 2. Example 1. Noisy case with multiple measurements. 
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Table 1. Example 1. The reconstruction errors 

11log /jl — log/i|| 2 / 11log/ j|| 2 for the noisy data. 


N 

1 

2 

3 

4 

5 

No regularization 

74.4% 

32.2% 

28.4% 

15.8% 

7.5% 

TV regularization 

23.2% 

6.0% 

22.1% 

4.8% 

4.6% 


It is evident that the larger N is, the more accurate the reconstruction becomes. 
More precisely, when N = 1 the reconstructed values of the absorption in the 
inclusions are completely wrong. This is due to the fact that the inclusions are 
roughly approximated by smooth atoms in A g and then corrected with fewer atoms 
in Af, and so the sparsest approximation does not separate the two components 
as desired. Choosing the same sparse approximations discussed above (so that 
CS2 is satisfied with N = 1) gives a big value for 5 in the noisy case, and so 
the reconstruction is not accurate. However, the problem is solved when more 
measurements are added: CS2 is easily satisfied with lower values of pf and p g 
when N becomes bigger (see Remark [l]). 

The reconstruction with N = 5 is very satisfactory if measurement and recon¬ 
struction errors are compared. Indeed, the noise from the data has almost disap¬ 
peared in the reconstruction, without a separate regularization. This is due to the 
implicit regularizing effect that sparse representations provide. 

We have also investigate the effect of an a priori total variation (TV) regulariza¬ 
tion [66] of the measurements hi on the reconstruction, using a Matlab implemen¬ 
tation based on the algorithm developed in jB4j • The regularization parameter was 
chosen a posteriori to achieve the best results, but in principle it can be learned 
from a training set [28] . The corresponding reconstruction errors with different 
values of N are shown in Table []] the improvement is significant only for a low 
number of measurements. For comparison, the regularized value of H<i is shown in 
Figure |2g| (where the usual staircase effect can be observed) and the reconstruction 
with N = 2 measurements is shown in Figure |2h[ 


5.2.2. Example 2 - The Shepp-Logan phantom. Here, we let ft be the well-known 
Shepp-Logan phantom (shown in Figure |3a|). We choose to stop the iterative proce¬ 
dure of OMP after 2000 iterations. As above, we consider the boundary conditions 
cpi as in (22) and the corresponding solutions Ui to (20), for i = 1,...,5, and 
measure the quantities hi as in (21). 


Table 2. Example 2. Reconstruction errors for the noise-free case. 


N 

1 2 3 4 5 

logM — log A 2 / 1 log A 2 

68.6% 24.8% 18.6% 11.4% 5.4% 


In a first experiment we consider the case without noise (Figure [3]). The recon¬ 
struction errors for N = 1,..., 5 are shown in Table [2] We see that the recon¬ 
struction quality improves as N increases, as it is expected from the general theory 
discussed in Section [3] From a comparison with the previous case without noise, 
we notice that more measurements are needed to have a satisfactory reconstruc¬ 
tion. This is due to the more complicated structure of the phantom, which has a 
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Figure 3. Example 2. Noise-free case with multiple measurements. 




HI-6 

1.4 
\ " 1.2 

I: 

(a) H 3 


(b) H 5 



(c) N = 5 



Figure 4. Example 2. Noisy case with N = 5 measurements. 


less sparse representation in terms of Haar wavelets than the absorption consid¬ 
ered in the previous case. Thus, a higher N is needed to satisfy the conditions in 
Definition |TJ 

In a second experiment we add white Gaussian noise ru to the data in (21). The 
noise level is such that 

, n ~ 17.8%. 

Motivated by the noisy-free case, we perform the reconstruction with N = 5 mea¬ 
surements. The reconstruction error is ||log fi — logjl\\ 2 / ||log/i|| 2 = 17.0%, that is 
comparable to the measurement error. The results are shown in Figure [4] It is 
expected that adding more measurements would improve the quality of the recon¬ 
struction. 
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5.3. Quantitative photo-acoustic tomography in the diffusive regime with 
variable T. 


5.3.1. Introduction. We consider here the problem of quantitative photoacoustic 
tomography, as introduced above, without the further assumption T = 1. Thus, 
the unknown absorption /i has to be reconstructed from 

H{pc) = T(x)ii(x)u{x), x E D. 

We consider the diffusion approximation ( p~9] ) of light propagation: 

(23) — div(Z)VR) + fm m 0 in Q. 


For simplicity, here we shall augment this equation with Dirichlet boundary condi¬ 
tions, that are supposed to be measurable. Note that D, T and /i are unknowns of 
the problem. In contrast with the case T = 1, where (23) was merely used to com¬ 


pute the data but not in the inverse problem, we shall make full use of this PDE. 
Let us briefly review the main known results on this inverse problem. Bal and Ren 
m showed that when T, D and /i are all unknown, then there is no uniqueness for 
the inverse problem even with all the measurements H for all solutions u to (23); 


namely, the PDE approach alone is not sufficient to reconstruct all the parameters. 
If any of these parameters is known, then the others may be reconstructed by using 
the PDE. The same authors have proved that all the coefficients may be uniquely 
reconstructed by using multi-frequency measurements, under certain assumptions 
on the dependency of the parameters on the frequency [2DJ. Naetar and Scherzer 
m studied the case of piecewise constant parameters: all the unknowns can be 
uniquely determined, but the method may be very sensitive to noise. 

We propose here for the single-frequency case a mixed approach combining the 
following aspects. 


As in m, the PDE ( [23] ) can be used in the reconstruction. However, one 
degree of freedom for the parameters remain. 

As in m, the PDE method gives unique reconstruction under the finite di¬ 
mensionality assumption of the coefficient spaces. 

The disjoint sparsity signal separation method may be applied to this case as 
in S[5T2] 


The combination of such approaches consists in substituting the piecewise constant 


assumption with the sparsity assumption, and then in the use of (23) to reconstruct 
all the parameters. More precisely, the reconstruction algorithm proposed here is 
substantially divided into the following three main steps. 

(1) By using the disjoint sparsity signal separation method applied to 

hi = log Hi = log(I» + loguj, i = 1,..., N, 

the solutions iq are reconstructed. 

(2) Following [T9] , by using the PDE 

u ■ 

—div(ThqV —) =0 in D, 

Ui 

with three suitable measurements, the diffusion D can be uniquely determined. 

(3) Finally, the absorption can be directly reconstructed via 

div(DViq) 


M = 


a 
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and possibly averaging over i. 


5.3.2. The reconstruction algorithm. Even though theoretically satisfactory, the al¬ 
gorithm summarized above is not applicable in practice as it stands. Indeed, the 
reconstruction of D in (2) is not too sensitive to errors in xxj, but that of /a in (3) 
is. To understand this, we compare the solutions xq to (23) with D = 1 and fi as in 
Figure 


la 


conditions given by (22): 


and the solutions u® to (23) with D = 1 and /x — /x° = 1, with boundary 


(24) 


—div(DVxq) + /xxq = 0 in $2, 

— div(D\/u^) + /a°u^ = 0 in fi, 
on dQ. 


Ui = U® = (fi 


The solutions are shown in Figure [5] Looking at the first two columns, it is clear 
that the variations between Ui and u® are minimal. This is due to the fact that 
the two problems have the same diffusion coefficients and small variations in the 
absorption coefficients. As we saw in § 5.2, the reconstruction at step (1) cannot 


be at this level of precision, and therefore fi cannot be reconstructed in this simple 
way. In order to overcome this difficulty, we make the following observation. 


Remark 4. The ratios xq/xx^ are almost independent of i, provided that /x is a small 
variation around a known background /x°. This is evident from the third column of 
Figure [5| and can be proved by arguing as follows. A direct calculation gives that 
Vi = Uiju\ satisfies 

f -div(WM + (M-M°M = 2^- mO, 

[ Vi = 1 on dfl. 

When fi is close to /x°, the right-hand side of this equation becomes negligible with 
respect to the other terms. Thus, Vi is substantially independent of i. 


Let us now describe the precise reconstruction algorithm based on these observa¬ 
tions. It consists of two initial steps and an iterative procedure consisting of three 
more substeps. For simplicity, we shall discuss only the noise-free case. We suppose 
that /x is a small perturbation around a known coefficient /x° and that D is known 
at one point of the domain. 

(1) By using the disjoint sparsity signal separation method applied to 

hi = log (T/x) + log Ui , i = 1,..., TV, 

a first approximation xq(0) of the solutions xq is reconstructed. As discussed 
in Section [4j this can be done by minimizing 

N 

mi n \\y\\ 0 +\Y^\\Afy f + Agyl-hiW 

ye R m f+ Nm 9 u y 

with OMP, and then writing xq(0) = exp (A g y l g ). 

(2) By using the computed xq(0) and the PDE 

77 * 

(25) — di v(Dui\7 — ) =0 in Q 

Ui 

with three suitable measurements, a first approximation T)(0) of the diffusion 
can be obtained. Indeed, choose three boundary values ipi, (f 2 and such 
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(a) ui 




(c) ui/u^ 



(d) U2 




0.994 


0.992 


(f) u 2 /ul 



(G) U 3 


Figure 5. 




(h) u° 

The solutions to 


24 


(0 U3/U3 

and their ratios. 


that 


(26) det(V—, V—) > 0, in n. 

Ui U\ 

(This can be easily done in two dimensions [19].) Then the above PDE may 
be rewritten as 


t (VlogD) = -[div( Wl V^)/ W2 div( Ml V^)Mi] [V£ v^] 


-1 


in 


which can be integrated over Q and gives a unique solution for the diffusion 
coefficient, since D is known at one point of the domain. Since the solutions 
Ui are very sensitive to changes in D , we expect this reconstruction to be 
satisfactory. From the numerical point of view, an optimal control approach 


may be applied to (25) to find D 


(3) We now start the main iterative procedure. Initialize /i(0) = and let rq(0) 
and T>(0) be as reconstructed in points (1) and (2). From /r(&), Ui(k ) and 
D(k), the following iteration is computed as follows. 
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(a) Given D(k) and /i(fc), let u^(k) be the unique solution to 

f -div(D(k)\7u o i (k)) + fi(k)u° i (k) = 0 inQ, 

\ u®(k) = tpi on dVt. 

Since ipi is known, u^(k) is a known datum. Therefore we can measure 
h°i = log (Hi/u^k)) = log(I» + log i = 1,... ,N. 

In view of Remark [ij the quantities Ui/u®(k) are almost independent of i. 
This leads to the minimization of 


N 

\\v\\o + *iY,W A fVf + A sVg- h i 

2/ GR m / + ( JV+1 ) m fl 




N 


X ^Y^\\ A f y f +A 9 y 9 +1 - h °i\ 


i= 1 


with Ai <C A 2 . The second term maintains the incoherence among the y l g s, 
on which this disjoint sparsity approach is based. The third term forces the 
quantities Ui/u^(k) to be independent of i, and numerical evidence shows 
that this gives a much better reconstruction than the one performed at 
point (1). The multipliers Ai and A 2 may be taken dependent on k. Set 
Ui(k+ 1 ) = u^(k) exp(A g y^ +1 ). 

(b) Given Ui(k + 1), find a better approximation D(k + 1) of the diffusion 
coefficient by proceeding as in (2). 

(c) Reconstruct fi(k + 1) via 

div(D(fc + +1)) infi 
v J N ^ Ui(k + 1) 


From the numerical point of view, it may be useful to regularize Ui(k + 1) 
and D(k + 1) before taking the derivatives. Finally, a TV-regularization 
of fi{k + 1) may reduce the accumulated noise. 

There is no obvious stopping criterion for this iterative procedure. However, in the 
numerical simulations less than five iterations were sufficient. 

In the above algorithm we have assumed for simplicity that the boundary values 
ifi are measurable. However, this is probably not a necessary conditions, since they 
may be obtained from point (1) as ipi = ^(O)^- 


5.3.3. Numerical simulations. We have tested the algorithm discussed above with 
the absorption map (1 considered in 


5.2. 1| (see Figure [6c]) . The same dictionaries 


considered in § |5.2| are chosen. The light intensities tq are the solutions of 


—di v(D\7ui) + jaui = 0 in £2, 
Hi = (fi on <9f2, 
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(d) D (e) /I (f) /z(2) 

Figure 6. The case with variable T. 


where the diffusion coefficient D is shown in Figure |6a| and five boundary values 
are chosen as follows: 

<Pl(x) = 1, 

(^ 2 (x) = 1 — sin(27rxi)/8, 

(^ 3 (x) = 1 - sin(27nr 2 )/8, 

(^ 4 (^) = aq/4 + 7/8, 
ip${x) = x 2 /4 + 7/8. 

The internal data take the form 

Hi = T/Hq, z = 1,... 5, 

where the Griineisen parameter is shown in Figure |6b| The measurements corre¬ 
sponding to the first three boundary conditions will be used for the disjoint sparsity 
signal separation method (with N = 3), namely for the steps (1) and (3a) of the 
above algorithm; those corresponding to (/?i, </? 4 and will be used in the steps (2) 
and (3b), in order to satisfy ( [26] ) . All measurements are used in the last step (3c). 
The OMP iterative procedures are stopped after 2000 iterations. If the absorp¬ 


tion /I were recovered via (27) immediately after steps (1) and (2), the reconstruc¬ 
tion would not be satisfactory, as it can be seen in Figure |6e[ This makes steps (3a) 
and (3b) necessary: after repeating step (3) twice, the quality of the reconstruc¬ 
tion is sensibly improved (see Figure [6f]). The corresponding reconstruction of D is 
shown in Figure [6d] As anticipated before, the reconstruction of D from rq is much 
more stable than that of fi. Note that the absorption was supposed to be known 
near the boundary of the domain, to avoid problems with the second derivatives in 


( 27 ) 
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The case with noise was not studied in this paper, since the robustness to noise 
of the disjoint sparsity signal separation algorithm has already been tested in § |5.2| 
The robustness to noise of the other steps of the reconstruction method discussed 
above is well-known, and standard regularization method can be employed. It is 
worth noticing that since the absorption g is found in step (3c) from the recon¬ 
structed values of the light intensities iq, the signal to noise ratio of iq has to be 
sufficiently large. Unfortunately, as it can be seen from the third column of Fig¬ 
ure |5j the amplitude of the information about g captured in iq is very small, and 
so the noise has to be comparable to this. 

6. Conclusions 

In this work we have studied a method for signal separation based on the disjoint 
sparsity of multiple measurements. A theorem giving unique and stable reconstruc¬ 
tion was proved. The result is based on the incoherence of the measurements. Then, 
the method was applied to hybrid imaging problems, and in particular to quanti¬ 
tative photoacoustic tomography. This technique has been successfully tested on 
several numerical simulations, and results to be very robust to noise. 

The incoherence between the measurements gi is the main foundation of the 
method, and the numerical simulations have shown that such property holds true 
with different solutions to the same PDE, which is the relevant case for hybrid 
imaging. It would be very interesting to prove this result in general. Randomly 
chosen boundary conditions may give the necessary incoherence. 

Robust Principal Component Analysis (rPCA) [31 j might be used as an alterna¬ 
tive to the disjoint separation method for the recovery of / and the g^s from the 
knowledge of hi ~ f + &, i = 1,..., N. Writing / = A f y f and gi = A g y l g , the 
known matrix M := t A g [hi • • • can be expressed as 

M = [AAA/ • • • AAA/] + [vl ■■■ Vg ], 

that is the sum of a rank-one matrix and of a sparse matrix. It would be inter¬ 
esting to see whether the requirements for the application of rPCA are fulfilled. 
Heuristically, this is indeed the case: the incoherence between Af and A g should 
provide the non-sparsity of the low-rank component, while the incoherence, or pos¬ 
sibly the disjoint sparsity, of the gi s should ensure that the sparse component is 
not low-rank. 

Finally, it would be very interesting to investigate whether the main ideas behind 
this method can be applied to other inverse problems with multiple measurements 
consisting of two components, of which only one remains fixed. 

Appendix A. Uncertainty principles 

A.l. Proof of Proposition [2j Let Af and A g be as in Proposition [2] For p = 
1,..., n and q G M n define 

£p(<?)= ^ max „ mm l<?A)l- 

Note that £ p (q) > 0 if and only if ||g|| 0 > p: in this sense, the map £ is a quantitative 
version of the norm ||-|| 0 . Let ceil(z) denote the nearest integer greater than or equal 
to z. 


DISJOINT SPARSITY FOR HYBRID IMAGING 


27 


Lemma 2. Take ( = 1,..., m g and p = 1,... , n such that 

11‘^vilo >P. veR m °\{0}, IMIo <c 

There exists C^> 0 depending only on n, Af and A g such that 

£p( AfA g v) > C^, 

for every v G st/c/i that \\v\\ 2 = 1 and ||p|| 0 < £. 

Proof. By contradiction, there exist vi G M 7 ™ 9 such that ||p/|| 2 = 1, ||pz|| 0 — C and 
0 as Z -A oo, where qi = t AfA g vi. Up to a subsequence, we have vi -A v 
for some v G M™ 9 such that \\v\\ 2 = 1 and ||p|| 0 < £ and qi -A l AjA g v =: g. By 
assumption we have ||g|| 0 > p. Thus, there exist 1 < oq < • • • < a p < n and 5 > 0 
such that |g(ag)| > 2s for every i = 1 , ...,p. Since qi -A g, there exists Zq such 
that |g/(c»^)| > 5 for every i = 1,... ,p and every Z > Zo- Hence f p (qi) > £ for every 
Z > Z 0 , a contradiction. □ 

We are now in a position to prove Proposition [2] 

Proof of Proposition^ Define D = ( v / m^+ IXl + C^ 1 ) and take g G M n such that 
||Ajg|| 2 > D and ||M.j-Hjg|| 2 < |. Write l A g Afq = v A r where 

v (a) = f A s A f q ^ if \ tA 9 A f<l( a )\ > r 
1 0 otherwise. 

Thus ||7*11^ < 1, whence ||r|| 2 < v Arif. Therefore, by the estimate 

D < \\ A fQ ll 2 = (IITA^IIa + - WAMU +1> 

and by construction of D we obtain 

(28) |M| 2 > - ||r|| 2 > D - Trig + jj) = (y/nig + |)C/ _1 > 0. 

Set v' := f/||p|| 2 , C := ||?/|| 0 and p = ceil(2/M) — (. By Lemma [ 2 ] whose 
assumptions are satisfied by Proposition [l] (using that Af l Af = /), we have 
tp^AfAgv') > C%. Thus, there exist 1 < aq < • • • < a p < n such that 

(29) | t AfA g v'(ai)\ > C^, for every i = 1.... ,p. 

From the identity Afq — Ag^AgAfq) + A^-^A^Afq), setting 2 = t Af(A g r A 
Ajj-^Ajj-Afq)) we obtain 

Q = C A f A g v ') IMI 2 + *■ 

In view of Lemma [l] parts 1 and 2 we have 

IML < Ml 2 ^ ||A r + A jC A j A fQ)\\ 2 < M 2 + < V™h + \- 

As a result, by ( [29] ) and ( [28] ) and the expression for D we have 

l?(“i)I = I IMI 2 C tA fA g v')(ai ) + z(ai)\ > |)«|| 2 Q - + |) > 0 

for every i = 1,..., p. Therefore 

Ikllo + IMIo > P + C > 2/M. 

Finally, the conclusion follows form the equality #{a : lO'AyAf(j)(o)\ > 1} = 
IRIlo- □ 
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A. 2. Proof of Proposition [3j Let us first discuss the orthobasis of 2D Haar 
wavelets of where d = 2 . Af is constructed via translations and dilations of 

four types of wavelets, as we now describe. Let j = 1 ,..., J — 1 denote the scale, 
from the finest to the coarsest and let £q, = 1 ,..., 2 J-J be the horizontal and 

vertical translation parameters. We consider four families of atoms ^ k defined by 


f ipj fe(2 J '(fc 1 - 1) + ai,2*(fc 2 - 1) + a 2 ) = -2-J, 

\ - 1) + o:i, 2 J (fc 2 - 1) + 2J- 1 + a 2 ) = 2"*', 

f ^ 2 fe (2 j (fci - 1) + a 1 ,2*(fc 2 - 1) + a 2 ) = -2"*, 

\ - 1) + 2J- 1 +a 1 ,2J(fc 2 - 1) + a 2 ) = 2"*, 

' - l) + 2i~ 1 +a 1 ,2i{k 2 - 1) + a 2 ) = -2“>, 

ipj k (2i(ki - 1) + ai, 2 J (fc 2 - 1) + 2J” 1 + a 2 ) = -2~\ 

' y 3 fe (2 J Ai - 1) + ai, 2 j (fe - 1) + a 2 ) = 2~1, 

k y 3 ’ fe (2^(fei - l)+2J- 1 +a 1 ,2J(fc 2 - 1)+ 2 J_1 +a 2 ) = 2"*, 

^ jfe (2- 7 (fci - 1) + ai,2- 7 (fc 2 - 1) +a 2 ) = 2 _J , 


Qq 1,... 

,2 J , 

G2 = 1, • • • 

12 J_1 , 

Gi = 1, . . . 

12 J '~ i , 

OL2 = 1, • • • 

,2 J ', 

Gi = 1, . . . 

to 

5-0. 

1 

G 2 = 1,... 

to 

t-o. 

1 

Qq , Qq = 1, 

• • • i 2 J 


and zero elsewhere. The orthonormal basis of Haar wavelets Af is given by 


J—2 

(J Mfc : i = 1,2,3, k e (1,..., 2 J “T 2 } U {VY I,k : * * 1, • • •, 4, fc e (1,2} 2 }. 

3 = 1 

The proof of Proposition [3] is based on the following result. 

Lemma 3. Under the assumptions of Proposition^ for every v G \ {0} we 
have 

J-B-l 

||*4/4 fl i,|| 0 > £ (2 J ~3 — 2L) 2 

3 = 1 


Proof Write g = By construction of A p , the vector g can be written as a 

linear combination of low frequency real sinusoids, namely 

^ = E E ytxi 

i= 1 h,l 2 = 0 

for some weights yf G M. Using standard trigonometric equalities, we can write the 
above sum in terms of complex sinusoids 

L 

5 (a) = E 0ie 2nih ^ e 2nil2 ^, a e { 1 ,..., 2 J } 2 . 

l u h=-L 

for some complex weights 61 G C. Since the constant vector is not in A g and v ^ 0, 
we have that 0i 7 ^ 0 for some l 7 ^ (0, 0). Without loss of generality, we can assume 
that there exists /* G {—L,..., L } 2 such that 

( 30 ) 0z* / 0 , r 2 ^ 0 . 

By construction of Af we have that 

J-B-l J-B-l 

(31) || *4/5110 > E #6{l,-,2 J - j } 2 :(5,^)2/0} =: £ #7} 

i=i i=i 

We now want to find a lower bound for fj=Tj. 
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Fix j = 1,..., J — B — 1. For simplicity of notation, set e^.(-) = e 27Tlli 2 J and 
rafc. = 2 J (/g — 1). For k G {1,..., 2 J_J } 2 we have 

L 2 j 2 j ~ 1 

23 Mlkh = E EE Oi (e/j (m kl + a\)ei 2 (m k2 + 2° 1 + a 2 ) 

ll ,l2 = — L Oi 1=1 Q!2 = l 

- ejj (m fcl + ai)e ^ 2 (m fc2 + a 2 )) 

L 2 j 2 J - 1 

= E d i e h( m ki)ei 2 (m k2 ) E E e G ai ) e * 2 ( a 2 )( e i 2 ( 2i_1 ) - !) 

l\^2 = — T Q!i = l Q 2-1 

L 

= E C h(ki)e l2 (m k2 ) 

h=~L 

where we have set 

L 2 j ~ 1 2 J ' 

6 2 (A)~ E ^( e / 2 ( 2 J - 1 )-l)(E e ^M)(E 

l\ — — L Oi 2=1 a\ = 1 

Since L <2 B with B < J — 2 and j < J—B — 1, by using standard identities for geo- 
metric sums it is easy to show that (e*^- 1 )- 1 ) (E« 2 =i e h( a 2)) (E« 1= i e h («i)) 7 
0 for all Zi and all I 2 7 ^ 0 . As a result, in view of ( |30| ) we have that the polynomial 
in the complex variable z 

pis= E _ 1 )(E e ^( a2 ))(E e h( a i)) zh+L 

l\ = — L Ct 2 = l Oi\=l 

is non trivial. By the fundamental theorem of algebra, it has at most 2 L zeros. 
Therefore, writing z = e 27 ™ mfc i 2 J , we obtain that ifEj > 2 J ~ J — 2L, where Ej = 
{fci G {1,... ,2 J-J } : 7 ^ 0}. Take now k\ G Ej. Arguing in a similar way 

and writing z = e 27 rzmfc 2 2 we have that 

L 

e L (m k2 )2 j (g,ipl k ) 2 = Y Ch( k i)z h+L - 

h=~L 

Since £/*(fci) 7 ^ 0 , this polynomial in 2 : is non trivial, and so has at most 2 L zeros. 
In other words, for k\ G Ej we have that 

#{fc 2 e {1,..., 2 J ->} : ( 5 , ^ fc ) 2 7 0} > 2 J - J ' - 2L. 

Recalling that ffEj > 2 J ~i — 2 L, this implies that fj=Tj > (2 J_ - 7 — 2 L) 2 . Finally, 
the result immediately follows from ( |3l| ). □ 

We are now ready to prove Proposition [3] The proof follows the same argument 
used for Proposition [2] 

Proof of Proposition^ Define D = |(1 + C^ 1 ) (where C £ is given by Lemma |2j) 
and take q G M n with ||A/g|| 2 > D and || h4^-A/g|| 2 < 2/3. Lemma[l]part 1 and 
the identity Ajq = Ag^AgAfq) + A^^A^Ajq) yield 

D < \\A f q\\ 2 = (|| t A g Afq\\ 2 2 + || 2 < || ^A/g|| 2 + 2/3 


DISJOINT SPARSITY FOR HYBRID IMAGING 


30 


whence by construction of D we obtain 

(32) \\ t A g A f q\\ 2 >2Cp/3. 

Set v = L AgAfq/\\ t A g Afq\\ 2 . By Lemma [ 3 J the assumptions of Lemma [ 2 ] are sat¬ 
isfied with p = ^^T-, B_1 (2 J_ - 7 — 2 L) 2 and ( = m g ; as a result ^p^AfAgv) > 
C In other words, there exist a± < ••• < a p such that \( t AfA g v)(ai)\ > 
C%. Morevoer, by Lemma [l] parts 1 and 2 we have || l AfA^ l A^Afq^ 2 < 2/3, 
whence || ^fA^A^Afq]]^ < 2/3. As a consequence, by ( [32] ) and the identity 
q— || t A g Afq\\ 2 l Af A g v + t AfA^ l A^Afq we obtain 

\qM > II 'AgAfqW^AfAgviai)] - I *A f Af 'AfA^ou) | > - | = 0, 

for every i = 1,... ,p. In other words, \\q\\ Q > p, as desired. □ 
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